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Abstract 

The phase diagram of the attractive hard-core Yukawa fluid derived previ- 
ously [M. Robles and M. Lopez de Haro, J. Phys. Chem. C 111, 15957 
(2007)] is used to obtain the liquid-vapor coexistence curve of real water. 
To this end, the value of the inverse range parameter of the intermolecular 
potential in the Yukawa fluid is fixed so that the ratio of the density at the 
critical point to the liquid density at the triple point in this model coincides 
with the same ratio in water. Subsequently, a (relatively simple) nonlinear 
rescaling of the temperature is performed which allows one to obtain the full 
liquid vapor coexistence curve of real water in the temperature-density plane 
with good accuracy, except close to the triple point. Such rescaling may be 
physically interpreted in terms of an effective temperature-dependent at- 
tractive hard-core Yukawa interaction potential which in turn introduces 
an extra temperature dependence in the equation of state. With the addi- 
tion of a multiplicative factor to obtain from the model the critical pressure 
of real water, the corrected equation of state yields reasonably accurate 
isotherms in the liquid phase region except in the vicinity of the critical 
isotherm and in the vicinity of the triple point isotherm. The liquid-vapor 
coexistence curves in the pressure-temperature and pressure-density planes 
are also computed and a possible way to improve the quantitative agreement 
with the real data is pointed out. 
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1. Introduction 

Due to its ubiquitous presence and importance both in our everyday life 
and in nature, water has been the subject of numerous experimental and 
theoretical studies[l]. Yet, despite all of this work, the full understanding 
of its properties is still far from complete and new efforts are called for. 

Out of the many interesting and puzzling aspects of the properties of 
water, a particular challenging and still open problem is the determination 
of the full phase diagram of the system. In fact, a stringent test of many 
theoretical models is often their predicted phase diagram. In this paper, we 
will consider an aspect of this problem, namely the possibility of using a 
relatively simple model to obtain a sensible description of the liquid phase 
and the saturation curve of liquid water. Concerning the fluid region, the 
International Association for the Properties of Water and Steam from time 
to time issues documents on their recommended formulation for the ther- 
modynamic properties of ordinary water substance for general and scientific 
use, the most recent oneQ] dated in September 2009 and related to the so 
called lAPWS Formulation 1995 j^. This formulation consists of a funda- 
mental equation for the Helmholtz free energy of water which is obtained 
from a sophisticated fitting procedure of a wealth of experimental results. 
It is claimed to be valid in the fluid region up to temperatures of 1273 K 
and pressures up to 1000 MPa and to represent very well the most accu- 
rate experimental data in the single phase region and on the vapor-liquid 
phase boundary. On the computational side, a rather successful model is 
the TIP4P modelj3| or its modified version the TIP4P/2005 modelfl. This 
is a simple rigid nonpolarizable model with a Lennard- Jones site on the 
oxygen and partial positive charges on the hydrogen atoms. Further, the 
negative charge is placed along the bisector of the H-O-H bond and the 
parameters in the potential are adjusted to match the experimental density 
and the enthalpy of vaporization of liquid water. The TIP4P/2005 model in 
particular provides an excellent description of the coexisting densities along 
the orthobaric curve but it fails to yield the correct critical pressure. As a 
final piece of information required for the purposes of this work, mention 
should be made here of the analytic equation of state of liquid water intro- 
duced by Jeffery and Austin j^. Such an equation of state is based on the 
principle of corresponding states of Ihm, Song and Mason [tI and on the free 
energy of strong tetrahedral hydrogen bonds proposed by Poole et al^. It 
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is accurate and probably useful for engineering calculations but contains 
many adjustable parameters which is clearly an undesirable feature. 

Recently [9i], using thermodynamic perturbation theory and a hard-sphere 
(HS) fluid as the reference system, we have derived an approximate equa- 
tion of state for the attractive hard-core Yukawa (HCY) fluid and com- 
puted the corresponding phase diagram in the three different thermody- 
namic planes. This model fluid, which has received a lot of attention in 



different contexts, [10| is interesting because in certain well defined limits 
its interaction potential has proved to yield thermodynamically equivalent 
results to those of the Lennard- Jones potential, the Coulombic interaction, 
the HS potential, or the adhesive HS potential. The main asset of our 
derivation lies in the fact that the results are analytical and have a rather 
simple form. What we want to do here is to take advantage of the above 
mentioned properties and use these results to obtain a reasonable descrip- 
tion of liquid water (including the vapor-liquid coexistence curve) with a 
simple model. In this way we want to provide, albeit with limitations, a 
positive answer to the question raised in the provocative title ( "Can simple 
models describe the phase diagram of water?") of a recent paper by Vega 
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The paper is organized as follows. In order to make it self-contained, 
in the next section we provide the equation of state of the attractive HCY 
fluid. This is followed in Sect. |3]with the description of the procedure to 
map the results of the attractive HCY fluid to those of real water in the fluid 
phase. Section H] is devoted to the mapping of the liquid- vapor coexistence 
curve in the pressure-density and pressure-temperature planes. The paper 
is closed in Sect. |5]with further discussion and some concluding remarks. 



2. Approximate equation of state for the attractive HCY fluid. 

The attractive HCY fluid is a system whose molecules interact via the 
pair potential 

'^^^^ ^""^ - \ -e^ exp[-^(r - a) /a], r>a ' 

where a is the hard-core diameter, e is the depth of the potential well at 
r = cr, and z is the (dimensionless) inverse range parameter. 

The potential (pHCY {r) may be split in a form convenient for the use 
of the liquid state perturbation theory as 4>hcy (r) = 0o (r) + 0i (r), where 
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the reference potential 0o will be a HS pair potential (where the "effective" 
diameter d of the HS reference system will simply coincide with the hard- 
core diameter, a), and the perturbed part 0i is given by 



^'^"^^ I -e^exp[-z\r-a)/a], r>a 



Next we consider the usual perturbative expansion [12| for the (excess) 
Helmholtz free energy Af^^jy of this fluid. To first order in /3 = l/ksT, with 
T the absolute temperature and ks the Boltzmann constant, its approxi- 
mate expression reads (for details see Ref. (oj) 

A ex A ex 

^MCX^ = _ i2r^i3ee^G{z), (3) 



where 



is the (excess) free energy of the HS fluid, Z^siv) = 1 + 4:V9Hs{.cr^) is the 
compressibility factor of the HS fluid, rj = is the packing fraction (p 

being the density and the number of particles), gHs{o'~^) is the contact 
value of the HS radial distribution function and G{z) = G{t = z) with 
G(t) = Clrgnsif^)], denoting the Laplace transform operator. Here gusij) 
stands for the radial distribution function of the HS fluid. 

Within the rational function approximation method[l^, G{t) turns out 
to be be given by 

where $(t) is a rational function of the form 

. i + Srt + s^e + s^e + Stt' 
*(*) = tTlJTlJ^ 

and the coefficients Si and Li are the following functions of the packing 
fraction t] 

1 ri + 12riL2 + 2 - 247^^4 

- 2 2^1 ' 

L2 = -3{Zhs-1)S^, (8) 
3 -1 + 4^9-8^4 



1 -r/ + 8t]L2 + 1 - 2L2 - 24r7^4 
"2 27] + 1 ' 

1 2r]-r]^ + 12r]L2{r] - 1) - 1 - 72r]^Si 
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1 — rj 



367] {Zhs- 1/3) 



{27] +1)7] 
1 - 



^ ^ Z^g - 1/3 / Xhs ^ 
Zhs ~ ZpY \Xpy , 



1/2 



(10) 

(11) 

>(12) 



^2 = 
^3 = 

^4 = 
with, ZpY , 

that the only input required to completely determine A'f^^y is an expression 
for Zhs- Here we will consider the compressibility factor that follows from 
the popular Carnahan-Starling equation of state given by 14 



Xpy 



and Xhs = {■^jivZHs]) ^ Note then 
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I +7] + 7]'^ 
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(13) 



which is known to be rather accurate in the stable fluid region. This in turn 
implies that in this case 



XHsiv) 



1 + 4?7 + 4r/2 - 4:7]^ + 



(14) 



From Eqs. (j3]) - fll3p one may easily derive the explicit form of the 
equation of state of the attractive HCY fluid in terms of the dimensionless 



variables P 
namely 



Pa^/e (P being the pressure), p = pa^ and T = ksTje, 



P = pTZ 



HS 



2vrpV|,[pG(^)], 



,6 V -^dp—- ^''^ 

where it is understood that the packing fraction r] appearing in all the 
expressions related to G{z) must be replaced by |p. 

The (approximate) phase diagram of the attractive HCY fluid in the 
three thermodynamic planes (T-p, P-T and P-p) was obtained using Eq. 
( 1T5|) and compared to simulation data in Ref. j9|. Such a comparison indi- 
cated that the approximate equation of state was a reasonable compromise 
between accuracy and simplicity. With this in mind, in the next section we 
will try to use the same equation of state to obtain the phase diagram of 
water in the stable fluid region. 
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3. Scaling for the thermodynamics of water 



The length and energy units for the attractive HCY fluid are deflned 
arbitrarily by the parameters a and e in the interaction potential. The 
parameter z plays the role of an adjustable free parameter needed in any 
application of the model. To map the attractive HCY fluid to the real 
thermodynamic properties of water, we have taken the following steps: 

• To flx the value of z in the way specifled below. 

• To flnd a non-linear rescaling of the temperature to flt the saturation 
branch of the coexistence curve. 

• To compare the rescaled attractive HCY fluid isotherms with those of 
real water. 

• To adjust the value of the critical pressure of the attractive HCY fluid 
to match the one of real water. 

3.1. The choice of the inverse range parameter 

In our previous work,(9(] we discussed how to obtain, for a given the 
liquid-vapor coexistence curve for the attractive HCY fluid through the 
usual Maxwell construction in the T-p plane. Now we want to determine 
the value of z such that the coexistence curve derived for the attractive HCY 
fluid may be used for water. For this purpose it is convenient to introduce 
the reduced quantities T* = T/Tc and p* = p/pt, where Tc and pt are the 
critical temperature and liquid density at the triple point, respectively. In 
what follows we will add to the subscripts the labels w and Y to indicate 
real water and attractive HCY fluid, respectively. For water, we take T^c = 



647.096 K and p^t = 999.793 kg/m^[l7| It should be pointed out that 
for our approximate equation of state for the attractive HCY fluid, Eq. 
(|T5l) . the density at the triple point is independent of z and is given by 



PYt = 0.9375. [15| On the other hand, the critical density and temperature 
for the attractive HCY fluid only depend on z. We next flx z by requiring 
that the ratio pyc/ PYt of the attractive HCY fluid coincides with the reduced 
critical density of water, pwc/ Pwt = 0.32207. In this way, the resulting value 
of the inverse range parameter is ^ = 1.7118 which yields TVc = 1-35997 
and pyc = 0.3019. Since Tyc = ksTc/ec (where we have added a subscript 
c to e for reasons that will become clear later) and pYt = PwtO'^, we get 
ec — 6.57 X 10^^^ J and a ~ 3.0401 x 10^^° m. Therefore we now have all 
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Figure 1: Coexistence liquid- vapor curves for real water (continuous line) and the fitted 
attractive HCY fluid with z = f .7118 (red dashed line). 



the quantities required for the comparison of the two coexistence curves. In 
figure [1] we show the hquid- vapor coexistence curves in the T*-p* plane for 
both real water and the corresponding attractive HCY fluid with the above 
value of z. 

Note that the shape of the two curves is similar. Nevertheless, one still 
has to manipulate the data if a better quantitative agreement is desired. 



Undoubtedly, the best way would be to follow a crossover treatment [18 



in- 



corporating the scaling laws valid for the critical region, as has been done for 



instance in the case of many compounds including water by Llovel et a/ [19 
in connection with a generalized van der Waals-type equation of state. We, 
however, have taken a much simpler (pragmatic) approach which suffices 
for our purposes and is described in what follows. 

3.2. Nonlinear rescaling for the temperature 

From the saturation branch of the coexistence curves in figure [1], we 
extracted the two different reduced temperatures corresponding to a fixed 
density and produced a plot T* vs Ty. This yields 



T* = f(t:). 



(16) 




Figure 2: Reduced temperature of the attractive HCY fluid as a function of the reduced 
temperature of water in the saturation branch of the vapor-liquid coexistence curve. 
The red dashed line corresponds to data extracted from figure [T] and the continuous line 
corresponds to Eq. (fT7|) . 



where a good fit of tlie data produced in tliis way is a non linear function 
of tfie form 

F{T:) = 1 - [ai(l - t:)*i/4 + a^il - T:)V3 + 

a,{l-T:y/' + a,{l-T:)], (17) 

with ai = 1.39076, oa = -3.34647, ag = 3.22505 and = -0.610568. 
Figure [2] displays the data corresponding to this fit. 

Now we address the question of whether the above information, namely 
Eq. fll7p . allows us to 'predict' the coexistence curve of real water in real 
units using the results of the HCY fluid shown in figure [H The answer 
proceeds as follows. For any temperature on the real water vapor-liquid 
equilibrium curve we first determine the corresponding and then use 
Eq. f[T7|) to obtain the corresponding Ty. Next we multiply it by Tyc and 
derive the corresponding vaporization and condensation densities pwap and 

PYcond- 

Finally the real densities are simply PwtPvvap/ Pvt and PwtPYcond/ Pvt- 
Schematically if we denote by fiTy, z) the mapping which allows one to 
determine {pyvap, PYcond) in our approximation for the attractive HCY fluid. 
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Figure 3: Liquid- vapor coexistence curve as obtained with the recipe ((H]) (dashed Une), 
compared to experimental data[l3 (continuous line) and the results of the TIP4P/2005 
model[l6[ (crosses). 




then the recipe to obtain the real water liquid-vapor equilibrium curve is: 

{pwvap, pwcond) = ^^/(-^(^to/Tii,c)7Vc, 1 .71 18) . (18) 
PYt 

Figure [3] shows the predicted coexistence curve using this recipe and the 
real water coexistence curve. In this figure we have also included, for com- 



parison, the results obtained from the TIP4P/2005 model. [16 

One can readily see the very good agreement between our prediction and 
the real coexistence curve, except close to the triple point. A particulary 
noteworthy feature is that Eq. f|T71) was derived by considering only the 
saturation branch, but it does a very good job also for the rest of the coex- 
istence curve. Moreover, Eq. f lT7|) admits the following appealing physical 
interpretation. Assume that, instead of considering a constant depth of the 
potential well at r = a (as we have done so far with e = ec), in order to 
make the HCY fluid thermodynamically equivalent to real water for each 
temperature T we are taking an effective HCY potential with the same a 
and z but with e = e{T). Then, Eq. f|T6|) can be rewritten as 

^/<^^>-f(L). (19) 



or, equivalently, 



In figures m and O we show the effective (temperature dependent) potential 
as a function of distance and the temperature dependence of its depth at 
r = cr, respectively. Note that the overall shape of the potential is preserved, 
but the effect of the extra (non monotonous) temperature dependence of e 
on the coexistence curve is clearly important. Although such a task lies 
beyond the scope of the present paper, we are persuaded that it would be 
interesting to investigate whether this effective interaction may provide an 
adequate picture of the real interaction between water molecules. 




Figure 4: Effective (temperature dependent) attractive HCY potential as a function of 
distance for different temperatures: Continuous line: T = 273 K\ dashed line: 473 K\ 
dotted line: 647 K. 



The next question that comes to mind is whether the same scaling of 
the temperature and density will hold for the equation of state. This is 
examined in the next subsection. 

3. 3. The equation of state of liquid water 

The explicit form for our approximation to the attractive HCY fluid 
equation of state is given in Eq. ( IT5l) but, with the above discussion, it will 
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Figure 5: Temperature dependence of the well depth e(T) as given by Eq. (PH)) . 



be understood that we shall use in it the temperature dependent e(T) as 
given by Eq. IHM . However, one finds that the critical pressure obtained 
in this way with the attractive HCY fiuid is Pyc = 38.08 MPa whereas the 
real critical pressure is P^c = 22.064 MPa. This overestimation seems to 
be due to a similar feature observed for the results derived within the ther- 
modynamic perturbation theory for the attractive HCY fiuid. In fact, in 
reference [oj we pointed out that the predicted critical pressure for the at- 
tractive HCY fiuid was always greater than simulation values. For instance, 
ii z = 1.8 (which is the value of the inverse range parameter closest to our 
z = 1.7118 here where simulation results are available), the critical pressure 
obtained from Eq. (|T5i) is P = 0.1576 whereas from simulations one gets 



approximately P = 0.101 |20|. This gives a ratio 0.101/0.1576 = 0.6409. In 
the present case, the ratio P^c/ Pyc yields 7 = 22.064/38.08 = 0.5794 which 
could be reasonably attributed to the discrepancy between the predicted 
critical pressure of the attractive HCY fiuid and the value that one would 
obtain from simulation z = 1.7118. Hence, we 'correct' our equation of 
state Eq. with a factor 7 = 22.064/38.08 = 0.5794 so that it now reads 
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6 op 



(21) 



Figure E] shows the result of using Eq. fl^ together with the experimen- 
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Figure 6: Isotherms as obtained from experimental data[13l (continuous lines), from Eq. 
(PT|) (dots) and from the Jeffery- Austin equation of state[6| (dashed hnes). From right 
to left the curves correspond to T = 300 K, T = 350 K, T = 400 K, T = 450 K, T = 500 
K, T = 550, K T = 600 K and T — Twc- The real coexistence line has been drawn to 
clearly indicate the liquid phase region. 



tal isotherms jl 71] and those that follow from the Jeffery- Austin equation of 
state j^. Clearly the agreement between this corrected equation of state 
and the experimental data for the isotherms is reasonably good in the liquid 
region except close to the critical and triple point isotherms. In the case of 
the critical isotherm, it leads to the correct critical point but deviates from 
the real critical isotherm at higher densities. On the other hand, the Jeffery- 
Austin results are also remarkably accurate, but show similar limitations. 
In particular they fail to produce the correct critical isotherm. 

4. Vapor-liquid equilibrium lines in the P — T and P — p planes. 

Once we have shown that the real vapor-liquid equilibrium coexistence 
curve may be accurately described in the T-p plane through the use of 
our approximate (corrected) equation of state, Eq. fl2T]) . it is reasonable 
to wonder if using the same approximate equation of state one may get 
also an accurate description in the P — T and P — p planes. The answer 
to the above question is only partially affirmative in the sense that the 
approximation is very rough. In order to get an accurate description one 
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Figure 7: Relationship between the reduced pressures of the real water and of the HCY 
fluid as derived from the coexistence curves in the P — T plane. The continuous line is 
the fit given in Eq. ^ 



is forced to make a subsequent rescaling on the pressure which, in this 
instance, is unfortunately not simple due to the form of the relationship 
between the reduced real pressure (P^ = Pw/Pwc) and the reduced HCY 
pressure {P* = Py/Pyc) derived from the coexistence curves in the P — T 
plane and given by the curve in figure [71 

In order to fit accurately this curve both in the low pressure end and 
close to the critical pressure it is unavoidable to introduce a high order 
polynomial. We have taken the following form: 

P'{P) = Y^c,P\ (22) 

i=l 

where the coefficients Cj are given by ci = —0.149421, C2 = 1.91855, C3 = 
-1.23117, C4 = 0.519336, C5 = -0.150666, Cg = 0.0306816, C7 = -0.004456, 
Cg = 0.000466, cg = -0.000035, Cw = 1.880416 x 10"^ Cn = -7.003285 x 
10-8, C12 = 1.719548 X 10-^ ci3 = -2.501064 x 10"" and Cm = 1.631257 x 
10-13. 

In figure [8], we show in a log-linear plot the comparison between the dif- 
ferent coexistence curves in the P-T plane. Note the disagreement between 
the experimental curve and the one obtained with the attractive HCY fiuid 
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Figure 8: Coexistence curves in the P -T plane in a log-linear plot. Continuous line: 



experimental data [171; dots: attractive HCY fluid with temperature rescaled using Eq. 
([T7|) : triangles: attractive HCY fluid with temperature rescaled using Eq. (pTl) and 
pressure rescaled using Eq. (22); crosses: results from the TPI4P/2005 modeljlq. 



using only Eq. fl2ip and the substantial improvement once the scaling im- 
plied by Eq. (12^ is introduced. The inclusion of the results obtained with 



the TPI4P/2005 model[l6[ serves to indicate that the qualitative tendency 
is reproduced by this model but the pressure is underestimated. 

Finally, we consider the P — p plane. Figure |9] contains the comparison 
between the different coexistence curves in this plane. Note once more the 
disagreement between the experimental curve and the one obtained with 
the attractive HCY fluid using Eq. (12T|) and the improvement when the 
scaling given by Eq. ^ is also used. Again, the TPI4P/2005 modelgQ 
correctly reproduces the qualitative trends but in particular the critical 
pressure is highly underestimated. 

5. Concluding remarks. 

As the results of the previous section indicate, we have been able to use 

a simple model of intermolecular interaction, the attractive HCY fluid, to 

obtain a very reasonable picture of the thermodynamic behavior of water 

in the stable liquid phase. This is remarkable given the fact that the model 

does not consider at all important aspects such as bond geometry or charge 
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Figure 9: Coexistence curves in the P - p plane. Continuous line: experimental data[l7[; 
dots: attractive HCY fluid with temperature rescaled using Eq. ([T7)) : triangles: attractive 
HCY fluid with temperature rescaled using E g. (TTTl) and pressure rescaled using Eq. ; 
crosses: results from the TPI4P/2005 model[l6|. 



distribution. In fact, our formulation involving an effective (temperature- 
dependent) attractive HCY interaction shares some aspects of the approach 
followed in colloidal systems, where some degrees of freedom are accounted 
for by an effective interaction potential. Moreover, the present develop- 
ment is yet another effort in devising a sound molecular theory of liquid 
water based on the principles of statistical mechanics and the thermody- 
namic perturbation theory of liquids. It is fair to say that with much less 
effort compared to the lAWPS Formulation 1995 |3i] and probably also to the 
Jeffery-Austin equation of state [6] we are able to achieve good accuracy. We 
should of course acknowledge the semiphenomenological character of our ap- 
proach that allowed us to have this description, namely the fit of the inverse 
range parameter of the model z and the manipulations involved in the vari- 
ous rescalings. Also worth mentioning in this regard is the limitation to the 
liquid phase, which becomes manifest in that the agreement between ex- 
perimental data and our formulation gets worse both near the critical point 
and near the triple point. Furthermore, the scaling implied in Eq. (|T7j) 
prevents us from dealing with the supercritical region which is of course in- 
teresting in many industrial processes. However, a similar procedure could 
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be devised if one wanted to deal with the equation of state of steam. In 
such a case, instead of rescahng the temperature as we have done here, 
one could rescale the density. This might also lead to a density-dependent 
pair potential which could be analyzed following the recent approach intro- 



duced by Zhou [2l|, although the use of density-dependent pair potentials 



is not free from controversy [22[. Such developments lie beyond the scope 
of the present paper. Nevertheless, we seem to have been able to capture 
the main physics behind the thermodynamic properties of liquid water over 
a wide range of both temperatures and pressures. Therefore, and given 
the fact that the value of z = 1.7118 that we used in our calculations is 
relatively close to the one {z ~ 2) in which the thermodynamic properties 
of the attractive HCY are similar to those of the Lennard- Jones fluid, we 
wonder whether the replacement of the Lennard- Jones centers by attractive 
HCY centers in models used in simulation such as the TPI4P/2005 modeljsf 
might improve the prediction say of the critical pressure. This may only be 
ascertained by performing the simulation and our hope is that, apart from 
the potential use by engineers of the above results, this work may find an 
echo in this latter respect. 
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